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PREFACE 


This  document,  funded  as  an  IDA  Central  Research  Project,  has  been  prepared  to 
examine  the  use  of  newly  developed  models  to  study  how  contamination  is  transported. 
Standard  use  of  streamlines  and  particle  tracking  methods  have  limitations.  However, 
recent  developments  in  numerical  simulations  of  oil,  water,  and  gas  flow  in  petroleum 

reservoirs  have  provided  the  opportunity  to  improve  the  speed  and  accuracy  of  current 
approaches. 

This  document  summarizes  research  on  streamlines  at  Stanford  University  and 
commercial  development  of  streamline-based  simulators  by  StreamSim  Technologies  for 
the  potential  aid  they  may  provide  in  environmental  cleanup.  Further  investigation  of 
such  methods  may  provide  enhanced  information  regarding  how  these  techniques  can 
best  be  applied  to  contaminant  transport  problems. 

This  document  was  researched  and  written  by  Martin  Blunt  and  Martha  Crane  of 
the  Department  of  Petroleum  Engineering  at  Stanford  University;  Rebecca  R.  Rubin  of 
the  Institute  for  Defense  Analyses  was  the  project  leader.  The  authors  wish  to  thank  Dr. 

Robert  L.  Hirsch  of  Advanced  Power  Technologies,  Inc.,  for  his  review  and  helpful 
suggestions. 
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THE  USE  OF  STREAMLINE-BASED  METHODS  TO  MODEL 
CONTAMINANT  TRANSPORT  IN  THE  SUBSURFACE 

A.  SUMMARY 

The  use  of  streamlines  lines  that  follow  the  instantaneous  flow  field — are 
standard  in  hydrology  to  map  capture  zones  of  wells  and  to  visualize  flow  patterns. 
Particle  tracking  methods,  where  contaminant  mass  is  transported  along  streamlines,  are 
also  widely  used.  These  methods  model  contaminant  transport  free  from  numerical 
artifacts,  but  their  application  is  limited  to  simple  single-phase  flow  problems.  When  the 
flow  field  changes  with  time,  when  non-aqueous  phase  liquids  and  air  are  also  present,  or 
where  there  are  complex  non-linear  geochemical  or  radioactive  interactions  between 
chemical  species,  particle  tracking  methods  have  severe  limitations.  For  such  cases, 
conventional  finite  difference  or  finite  element  numerical  models  are  used.  But  these 
suffer  from  numerical  errors  and  grid  orientation  effects  and  tend  to  be  veiy  slow. 

Recent  developments  in  the  numerical  simulation  of  the  flow  of  oil,  water,  and  gas 
in  petroleum  reservoirs  have  extended  the  use  of  streamline-based  methods  to  complex, 
non-linear  problems  that  are  mathematically  similar  to  those  encountered  in  contaminant 
transport.  This  offers  the  opportunity  to  develop  streamline-based  methods  in  hydrology 
and  environmental  engineering,  resulting  in  a  significant  improvement  in  speed  and 
accuracy  compared  with  current  simulation  approaches. 

In  this  report,  we  first  give  a  brief  summary  of  research  on  streamlines  at  Stanford 
University  and  the  commercial  development  of  streamline-based  simulators  by  StreamSim 
Technologies.  We  then  provide  a  review  of  streamline  methods  and  particle  tracking  in 
the  hydrology  literature,  describe  new  ideas  in  a  streamline-based  approach  to  simulating 
contaminant  transport,  and  outline  its  potential  applications. 
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B.  RESEARCH  AT  STANFORD  AND  COMMERCIAL  DEVELOPMENT 

Much  of  the  recent  research  on  streamline-based  simulation  for  predicting  oil 
recovery  in  petroleum  reservoirs  has  been  performed  at  Stanford  University  in  the 
Department  of  Petroleum  Engineering  by  Rod  Batycky,  Marco  Thiele,  and  Martin  Blunt. 
Current  research  at  Stanford,  performed  by  Martha  Crane,  is  exploring  the  use  of 
streamline-based  methods  to  model  contaminant  transport.  This  work  is  being  performed 
in  collaboration  with  Andy  Tompson  and  co-workers  at  Lawrence  Livermore  National 
Laboratory  with  the  aim  of  applying  the  method  to  study  the  movement  and  fate  of 
radionuclides  at  the  Nevada  Test  Site. 

In  the  oil  industry,  there  is  considerable  interest  in  commercial  applications  of 
streamline-based  reservoir  simulation.  StreamSim  Technologies  was  founded  in  1997  by 
Rod  Batycky,  Martin  Blunt,  and  Marco  Thiele  to  develop  a  streamline-based  reservoir 
simulator  based  on  research  conducted  at  Stanford  Urfiversity  from  1991  to  1996.  The 
development  effort  is  supported  by  a  consortium  of  major  oil  companies.  The  primary 
application  of  the  simulator  will  be  to  predict  oil  recovery  from  detailed  geological  models 
of  the  reservoir  using  significantly  less  computer  time  than  traditional  finite-difference 
approaches.  The  technology  is  particularly  suitable  to  analyze  “what  if’  scenarios,  such  as 
multiple  geological  models,  infill  drilling  benefits,  pattern  conversions,  and  enhanced 
recovery  options.  In  the  future,  we  anticipate  that  a  commercial  contaminant  transport 
simulator  could  be  developed  using  the  same  ideas  and  that  could  be  used  to  predict  the 
fate  and  movement  of  contaminants  in  the  subsurface  and  aid  in  the  design  of  optimal 
remediation  and  containment  strategies. 

1.  Introduction 

Modeling  of  contaminant  transport  has  been  the  focus  of  much  attention  over  the 
past  few  decades  because  the  ability  to  make  quantitative  predictions  about  flow  and 
transport  in  the  subsurface  is  needed  to  properly  utilize  ground-water  resources  and  to 
restore  polluted  ground  water.  The  ability  to  do  realistic  ground-water  flow  simulation 
requires  a  quantitative  description  of  the  hydrogeological  setting,  which  often  requires 
transforming  more  qualitative  geological  data  into  single  or  multiple  realizations  of  the 
hydraulic  properties.  Although  there  is  uncertainty  in  the  assignment  of  hydraulic 
properties  throughout  the  flow  field,  there  is  a  consensus  that  ground-water  flow  modeling 
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has  been  fairly  successful  at  determining  an  accurate  pressure  distribution  and  velocity 
field.  For  example,  models  which  solve  the  ground-water  flow  equation  numerically  using 
either  finite  difference  or  finite  element  techniques  have  yielded  good  results. 

When  the  migration  of  dissolved  constituents  by  ground-water  flow  is  considered, 
the  difficulty  of  simulation  is  greatly  increased  because  of  the  complexity  of  the  transport 
process.  The  transport  of  solutes  in  ground  water  is  effected  by  a  large  number  of 
processes,  including  advection,  dispersion,  diffusion,  chemical  processes,  and  biochemical 
processes.  Because  of  this  wide  range  of  chemical  and  physical  processes,  scientists  from 
many  disciplines  have  contributed  to  the  understanding  of  contaminant  transport  in  the 
subsurface.  In  fact,  solute  transport  modeling  is  now  recognized  as  an  interdisciplinary 
challenge  (Abriola,  1987).  The  mathematical  statement  of  the  transport  of  a  non-reactive 
solute  is  the  advection-dispersion  equation: 

£-V.(D.VC)-V.(vC)  (D 

where  C  is  solute  concentration,  v  is  the  Darcy  velocity  of  the  ground  water,  and  D  is  a 
dispersion  coefficient.  Many  numerical  methods  for  solving  the  advection-dispersion 
equation  have  been  explored  by  ground-water  hydrologists. 

This  report  will  provide  an  overview  of  the  methods  most  commonly  used  in 
modeling  the  transport  of  solutes  in  ground  water.  Such  models  consider  advective, 
dispersive,  and  reactive  transport.  Particular  attention  will  be  given  to  the  ground-water 
literature  concerning  advection-dominated  transport  because  it  is  generally  agreed  that  for 
heterogeneous  porous  media,  especially  if  there  are  strong  sources  and  sinks,  considering 
advective  forces  only  provides  an  adequate  approximation  (Schafer-Perini  and  Wilson, 
1991).  In  fact,  in  field  cases,  macroscopic  dispersion  is  believed  to  be  controlled  by  the 
heterogeneity  of  hydraulic  properties  (Smith  and  Schwartz,  1980).  Therefore,  if  the  spatial 
variability  of  properties  such  as  hydraulic  conductivity  is  sufficiently  described,  ignoring 
dispersive  transport  is  a  reasonable  approximation.  Finally,  streamline  methods  developed 
in  the  petroleum  literature  will  be  described  and  compared  with  the  particle  tracking 
approach  m  order  to  consider  future  directions  for  contaminant  transport  modeling. 

2.  Numerical  Methods  for  Modeling  Solute  Transport 

In  the  ground-water  literature,  transport  models  are  categorized  as  Eulerian, 
Lagrangian,  or  mixed  Eulerian-Lagrangian.  In  the  Eulerian  approach,  the  transport 
equation  is  solved  on  a  fixed  spatial  grid,  so  concentrations  are  associated  with  fixed 
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points  or  volume  elements  in  space  (Bear,  1972).  Both  finite  difference  and  finite  element 
methods  are  examples  of  this  method.  Because  of  their  wide  use  and  success  in  flow 
simulations,  these  methods  were  among  the  first  used  in  transport  modeling.  It  is  now 
recognized  that  these  methods  handle  dispersion-dominated  transport  accurately  and 
efficiently;  but,  in  advection-dominated  transport  problems,  Eulerian  methods  suffer  from 
numerical  dispersion  and  artificial  oscillations,  especially  in  the  region  of  sharp 
concentration  fronts  (Kinzelbach,  1986;  Bear  and  Verruijt,  1987).  In  order  to  minimize 
numerical  errors,  small  discretization  in  time  and  space  is  needed,  which  requires 
enormous  and  sometimes  prohibitive  computational  effort  especially  when  considering 
field  scale  problems.  Since  many  field  situations  are  advection  dominated,  an  alternative 
to  Eulerian  methods  was  needed  to  handle  field  scale  solute  migration  in  heterogeneous 
porous  media. 

In  Lagrangian  methods,  concentration  is  associated  with  fluid  elements  or  particles 
that  move  with  the  prevailing  velocity  field  (Bear  and  Verruijt,  1987;  Zheng  and  Bennett, 
1995).  This  method  avoids  directly  solving  the  advection-dispersion  equation  by 
representing  the  solute  mass  by  a  large  number  of  particles  which  move  with  the  ground- 
water  velocity.  The  Lagrangian  formulation  of  the  advective  transport  equation  is 
presented  in  detail  in  Zheng  and  Bennett  (1995),  but  the  final  result  can  be  given  as: 
DC/Dt—O,  where  C  is  the  concentration  of  a  particular  fluid  element.  This  equation 
reflects  the  notion  that,  as  a  fluid  particle  moves  along  its  pathline,  its  concentration, 
provided  the  transport  is  purely  advective,  does  not  change.  Therefore,  the  solution  to  the 
advective  transport  problem  is  reduced  to  defining  pathlines.  Lagrangian  methods  are  free 
of  numerical  dispersion  and  are  accurate  and  efficient  for  modeling  advection-dominated 
transport.  These  methods  have  been  quite  successful  in  representing  the  movement  of 
steep  concentration  fronts  (Moltyaner  et  al.,  1993).  Because  of  the  success  of  the 
Lagrangian  approach  for  solving  the  advective  transport  equation,  these  methods  are 

widely  used  for  modeling  field  scale  advective  transport  (e.g.,  Bair  et  al.,  1990;  Guven  et 
al.,  1992). 

When  dispersive  transport  must  be  considered,  mixed  Eulerian-Lagrangian 
methods  (such  as  the  method  of  characteristics  approach  in  Konikow  and  Bredehoff  s 
(1978)  widely  used  two-dimensional  solute  transport  model)  have  attempted  to  combine 
the  strengths  of  the  two  previous  methods.  The  Lagrangian  approach  is  used  to  solve  the 
advection  term  of  the  transport  equation,  and  the  Eulerian  approach  is  used  for  the 
dispersive  term.  In  formulating  a  solution,  first  the  advective  transport  is  solved  using  a 
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Lagrangian  approach  in  which  each  particle  is  assigned  an  initial  concentration  value. 
Particles  are  moved  along  pathlines  for  one  time  step  and  placed  in  a  gridblock  on  the 
underlying  Cartesian  grid.  After  an  advective  step,  each  cell  node  is  assigned  a 
concentration  by  averaging  the  concentration  of  all  particles  within  that  cell.  Changes  in 
concentration  due  to  dispersive  transport  are  then  calculated  using  a  finite  difference 
method  on  the  rectangular  grid.  After  this  step,  the  particle  concentration  values  are 
updated  to  reflect  the  changes  in  grid  concentration  values.  Then  another  advective  step  is 

taken.  Zheng  (1993)  extended  the  method  of  characteristics  to  model  three-dimensional 
solute  transport. 

Moltyaner  et  al.  (1993)  compared  the  different  methods  for  numerical  simulation 
of  tracer  transport  in  a  field  scale  experiment  at  the  Twin-Lake  site.  The  results  indicated 
that  while  the  finite  element  model  chosen  for  this  study  could  successfully  describe  the 
pressure  field,  the  transport  solution  suffered  from  numerical  dispersion.  For  both  the 
method  of  characteristics  approach  and  a  Lagrangian  method,  the  random  walk  particle 
method  (Prickett,  1981),  superior  results  in  the  transport  simulation  were  obtained.  The 
final  assessment  in  this  paper  was  that  the  best  results  were  achieved  by  the  random  walk 

particle  method;  and  the  subsequent  studies  of  tracer  flow  at  the  site  used  this  approach 
only. 


a.  Modeling  Advective  Transport 

In  many  field  cases,  or  as  a  first  approximation,  transport  is  considered  to  be 
dominated  by  advection  (Bear  and  Verruijt,  1987;  Zheng  and  Bennett,  1995).  For  this 
reason,  many  studies  of  solute  migration  have  focused  on  purely  advective  transport 
(Guven  et  al.,  1992;  Bair  et  al.,  1990).  When  dispersion  is  neglected,  ground-water  flow 
paths  coincide  with  the  paths  of  contaminant  solutes.  These  pathways  have  typically  been 
determined  using  particle  tracking.  An  alternative  to  particle  tracking  for  defining 
pathlines  is  stream  functions.  Stream  functions  have  been  used  in  two-dimensions  to 
assess  the  travel  time  of  contaminants  (Nelson,  1978;  Javadel  et  al.,  1984;  Fogg  and 
Senger,  1985;  Frind  and  Matanga,  1985).  Approximating  streamlines  in  three-dimensions 

by  using  stream  functions  has  also  been  attempted,  but  the  mathematics  are  difficult 
(Strack,  1984). 


series  of  time  steps.  A  concentration  value  is  assigned  to  each  fluid  particle. 
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contaminant  transport  modeling,  non-contaminant  particles  receive  a  concentration  of 
zero;  contaminant  particles  receive  a  finite  concentration  so  that  the  entire  solute  mass  is 
represented.  Clearly  the  mass  assigned  to  each  particle  is  determined  by  the  number  of 
particles  chosen  to  represent  the  solute  species.  The  contaminant  particles  are  distributed 
either  randomly  or  uniformly  in  the  source  region  (Zheng  and  Bennett,  1995).  The 
particles  are  then  advanced  through  the  region  of  interest  along  path  lines  calculated  from 
the  flow  field.  The  description  of  the  flow  field  can  be  analytical  in  which  case  the 
velocity  is  known  everywhere;  but,  more  commonly,  the  velocities  are  defined  on  a  grid 
used  to  solve  the  discretized  flow  equations.  When  the  values  of  velocity  are  known  only 
at  the  interfaces  of  gridblocks,  an  interpolation  scheme  is  necessary  to  evaluate  the 
velocity  at  any  point  in  the  domain. 

One  major  difference  in  the  particle  tracking  methods  described  in  the  ground- 
water  literature  involves  the  choice  of  interpolation  schemes.  The  most  commonly  used 
velocity  interpolation  schemes  for  particle  tracking  are  simple  linear  interpolation  or 
multilinear  interpolation.  In  linear  interpolation,  each  component  of  the  velocity  vector 
varies  linearly  in  its  own  coordinate  direction.  Therefore  changes  in  the  x-velocity 
component  are  independent  of  changes  in  the  y  and  z  directions  (Pollock,  1988). 
Multilinear  interpolation  considers  each  component  of  the  velocity  vector  to  be  a  linear 
function  of  all  the  coordinate  directions.  For  example,  in  bilinear  interpolation,  which  is 
used  for  two-dimensional  models,  the  x-component  of  the  velocity  is  formulated  as  a 
linear  function  of  both  the  x  and  y  positions  (Anderson  and  Woessner,  1992).  These  two 
interpolation  schemes  give  somewhat  different  representations  of  the  flow  field.  In 
essence,  linear  velocity  interpolation  satisfies  a  cell-by-cell  mass  balance  but  gives  a 
discontinuous  velocity  field,  and  the  multilinear  schemes  give  a  completely  continuous 
velocity  field  but  fails  to  satisfy  mass  balance.  Other  velocity  interpolation  schemes,  such 
as  bicubic  interpolation,  have  been  used  as  well  but  are  not  common  (Zheng  and  Bennett, 
1995).  The  choice  of  a  velocity  interpolation  scheme  hinges  on  the  numerical  method 
used  to  solve  the  flow  problem.  If  a  finite  difference  method  is  used  to  solve  for  the 
velocity  field,  only  a  simple  linear  interpolation  scheme  is  consistent  with  this  formulation 
and  will  therefore  conserve  mass  locally  within  each  finite  difference  cell  (Pollock,  1988). 
Pollock  s  method  has  been  extended  and  applied  to  finite  element  representation  of  the 
velocity  field  (Cordes  and  Kinzelbach,  1992). 
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Choosing  a  velocity  interpolation  method  is  the  first  step  in  particle  tracking.  The 
next  and  final  step  involves  choosing  a  tracking  scheme.  Particles  are  tracked  along 
pathlines  by  solving  the  equations: 


II 

1 

(2) 

1 

II 

(3) 

ii 

1 

(4) 

The  methods  that  are  typically  used  to  solve  these  equations  include 
semianalytical,  Euler,  and  Runge-Kutta  (Anderson,  1992).  The  semianalytical  solution  is 
outlined  in  Pollock  (1988)  and  is  possible  only  if  a  linear  velocity  interpolation  scheme  is 
used.  Euler’s  method  involves  the  simple  numerical  integration  (e.g.,  in  the  jc-direction): 

xn+J  =  xn  +  vjxn,y„,z„)dt.  (5) 

This  was  the  tracking  approach  used  in  some  of  the  early  particle  tracking  literature 
(Konikow  and  Bredehoft,  1978;  Prickett,  1981).  Since  the  value  of  velocity  at  the  starting 
point  is  extrapolated  over  the  entire  interval,  the  time  step  must  be  sufficiently  small  to 
make  accurate  predictions  about  the  particles’  movement;  for  this  reason,  it  may  take 
several  time  steps  to  advance  a  particle  through  one  gridblock.  Euler’s  method  is  only  a 
first  order  approximation  of  the  integral,  so  the  accuracy  has  been  improved  by  using 
higher  order  numerical  integration  methods  such  as  the  fourth-order  Runge-Kutta  method 
(Zheng,  1989).  Euler  s  and  Runge-Kutta  methods  can  be  used  with  any  velocity 
interpolation  scheme.  However,  since  a  linear  interpolation  scheme  is  consistent  with  a 
finite-difference  approximation  of  the  velocity  field,  and  the  semianalytical  tracking 
method  introduces  no  numerical  error.  Pollock’s  method  for  pathline  computation  is  the 
best  for  using  with  finite-difference  flow  models. 

b.  Semianalytical  Pathlines 

Because  Pollock  s  method  for  pathline  generation  is  most  commonly  used  because 
of  its  appropriateness  for  interpolating  a  finite-difference  generated  velocity  field,  the 
details  of  this  tracking  procedure  are  included.  As  mentioned  before,  the  pathline  each 
individual  fluid  particle  follows  is  determined  from  knowledge  of  the  velocity  field  and  an 
underlying  assumption  that  the  velocity  field  varies  linearly  in  each  coordinate  direction 
and  is  independent  of  the  velocities  in  the  other  directions.  The  velocity  in  the  x  direction 
is  defined  as 

Vx=VX0+mx(x-xo)  (6) 
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where  mx  is  the  velocity  gradient  defined  as 

V,.*  - 


m  —  ■ 


Ay 


(7) 


Knowing  that  v,  =  dx/dt.  Equation  6  can  be  integrated  to  find  the  time  required  to 
reach  the  x  exit  face  Atex, 


A  /  =  — 


Vx.o  +  mx(xe  -  x„) 
,[Vx.o  +fnAX,-Xtl) 


(8) 


The  time  to  exit  the  other  faces  is  derived  in  similar  fashion.  The  pathline  will  exit 
the  face  with  the  smallest  exit  time  value  (Ate).  Finally,  the  exit  position  is  computed  by 
substituting  Ate  into  Equation  8  and  solving  for  xe, 


(9) 


When  a  pathline  is  traced,  the  inlet  and  exit  position  in  each  gridblock  as  well  as 
the  time  to  cross  each  gridblock  (Atej)  is  recorded.  The  gridblock  k  containing  the  particle 
traveling  on  a  particular  pathline  at  time  t  is  then  calculated  by  summing  the  time  it  takes 
the  pathline  to  cross  each  gridblock  until  the  time  of  interest  is  reached: 

i~k  /=*+] 

XA^./ -r- (io) 


/=o 


/= o 


It  is  important  to  point  out  what  information  emerges  from  a  particle  tracking 
approach.  The  main  results  of  particle  tracking  methods  are  mass  arrival  time,  arrival 
position,  and  concentration.  Because  these  are  the  data  required  by  many  environmental 
agencies,  particle  tracking  is  commonly  used  to  model  regulatory  problems  involving 
ground-water  contamination  (Anderson,  1995). 

Particle  tracking  has  been  used  to  study  the  pattern  and  rate  of  ground-water 
movement  and  contaminant  transport.  In  particular,  it  has  been  used  to  delineate  regional 
flow  systems  as  well  as  recharge  and  discharge  area  for  the  ground-water  system  of  Long 
Island,  New  York  (Buxton  et  al„  1991).  The  results  from  this  study  indicated  that  particle 
tracking  yields  results  that  are  consistent  with  the  conceptual  understanding  of  the  system. 
Particle  tracking  has  also  been  used  to  delineate  the  capture  zone  of  wells  (Shafer,  1987; 
Bair  et  al.,  1990;  Schafer-Perim  and  Wilson,  1991).  In  general,  capture  zones  of  wells  are 
delineated  through  reverse  particle  tracking  where  particles  are  started  at  the  well  and 
traced  backward  to  a  source  region.  Schafer-Perini  and  Wilson  (1991)  also  suggested  a 
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method  for  dynamically  allocating  particles  in  order  to  best  capture  the  migration  of  the 
contaminant  front.  Bair  et  al.  (1990)  used  particle  tracking  to  predict  the  pathways  of 
contaminants  from  hypothetical  release  points  along  a  highway.  Guven  et  al.  (1992)  used 
particle  tracking  to  model  a  two-well  tracer  test.  Their  model  considered  only  advection 
when  describing  the  transport  of  tracer  between  an  injection  and  production  well.  The 
results  of  this  study  indicated  that  purely  advective  models  were  successful  when 
sufficient  knowledge  of  the  spatial  distribution  of  hydraulic  conductivity.  In  recent  years, 
some  of  the  more  interesting  applications  of  particle  tracking  involve  coupling  with  other 
models  to  explore  the  effects  of  such  mechanisms  as  dispersion  (Prickett,  1981)  and 
geochemistry  (Fabriol  et  al.,  1993). 

c.  Modeling  Other  Transport  Processes 

One  major  limitation  of  particle  tracking  is  that  it  simulates  only  advective 
transport.  Some  work  has  been  done  to  incorporate  particle  tracking  into  more  complete 
transport  models.  Most  importantly,  the  effects  of  dispersion  on  contaminant  transport 
were  introduced  in  the  “random  walk”  model  (Prickett'  1981).  The  random  walk  method 
uses  the  particle  tracking  technique  of  associating  a  mass  of  solute  with  each  particle,  and 
then  the  effect  of  dispersion  is  included  by  adding  a  random  displacement  to  the  particle 
location  after  each  advective  time  step.  Basically,  dispersion  can  be  incorporated  in  a 
particle  tracking  method  by  adding  a  random  motion  to  the  motion  along  pathlines.  For 
the  spreading  of  solute  by  dispersion,  the  random  dispersive  displacement  can  be 
described  by  a  normal  distribution.  Since  Prickett  first  introduced  the  use  of  random  walk 
to  model  the  advection-dispersion  equation,  numerous  other  researchers  have  used  the 
method  to  look  at  transport  in  porous  media  (e.g.,  Uffink,  1988;  Gelhar,  1990;  LaBolle  et 
al.,  1996).  Dispersion  is  an  important  mechanism  to  model  in  solute  transport  because  a 
dispersed  plume  is  more  widespread  than  a  plume  moving  by  advection  alone;  also  the 
concentration  of  solute  in  the  plume  is  reduced  by  the  dispersive  process  (Mercer  and 

Waddell,  1993).  In  addition,  the  arrival  of  the  first  contaminant  particles  at  a  point  of 
interest  is  often  controlled  by  dispersion. 

Simple  chemical  reactions  have  also  been  included  in  particle  tracking  codes.  For 
example,  adsorption  has  been  modeled  by  adding  a  retardation  factor  to  the  mass  transport 
equation  (Wen  and  Kung,  1995).  The  process  of  adsorption  reduces  the  rate  at  which  the 
contaminant  front  is  moving  because  of  reactions  between  the  solid  matrix  and  the  solute 
which  cause  some  of  the  solute  to  bond  to  the  solid  surface.  If  this  reaction  is  fast  (relative 
to  the  physical  transport),  reversible,  and  modeled  by  a  linear  isotherm,  the  effect  can  be 
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represented  by  a  retardation  factor.  The  retardation  factor  allows  the  transport  equation  to 
be  expressed  in  a  purely  advective  form  by  replacing  the  Darcy  velocity  with  a  retarded 
velocity  given  by: 

v 

^ retarded  D  ’  where  R  is  a  retardation  factor. 

R 

This  is  an  important  mechanism  to  understand  because  adsorption  makes  it  more 
difficult  to  remove  the  solute  at  a  given  site.  Another  process,  which  has  been  included  in 
particle  tracking  codes,  is  radioactive  decay  (Wen  and  Kung,  1995).  Obviously,  this  is  an 
important  phenomenon  to  consider  at  sites  polluted  with  nuclear  wastes.  If  the  half-life  is 
of  the  same  order  or  less  than  the  residence  time  in  the  ground-water  system,  then 
radioactive  decay  can  be  an  important  mechanism  for  attenuation.  Also,  the  daughter 
products  of  the  radioactive  species  add  to  the  complexity  of  the  system  (Mercer  and 
Waddell,  1993).  Radioactive  decay  is  included  in  particle  tracking  analysis  through  the 
addition  of  a  decay  constant  (1)  to  the  mass-transport  equation.  (Kinzelbach,  1986). 

Essentially,  the  mass  associated  with  each  particle  changes  over  time  according  to  the 
relation: 

^ panicle  ~  ^panicle  e  M' ,  where  M  is  the  mass  of  the  particle. 

Goode  and  Konikow  (1989)  modified  the  method  of  characteristics  code  to  account 
for  decay  as  well  as  equilibrium  controlled  sorption. 

d.  Assessment  of  the  Particle  Tracking  Approach 

Clearly,  the  Lagrangian  approach  has  many  advantages  for  modeling  contaminant 
transport.  The  problem  of  numerical  dispersion  is  eliminated.  Some  of  the  key  transport 
processes  can  be  included  in  a  straightforward  and  intuitive  way.  But,  there  are  drawbacks 
to  the  method  as  well.  The  concentration  field  is  computed  by  evaluating  the  location  and 
mass  of  the  particles.  The  accuracy  of  this  method  depends  completely  on  the  number  of 
particles  chosen  to  represent  the  contaminant  mass.  Representing  more  than  one  solute 
species  requires  a  huge  number  of  particles;  and  the  method  loses  its  computational 
efficiency.  Because  of  these  difficulties,  modeling  multi-species  transport  and  considering 
interactions  or  coupling  between  species  is  not  possible.  Since  most  contaminants  are 
composed  of  many  components  that  interact  with  each  other,  particle  tracking  is  somewhat 
limited  in  its  application  to  practical  problems. 

It  has  become  clear  that  for  solute  transport  modeling  to  be  a  predictive  tool,  the 
reactive  processes  that  affect  the  movement  of  contaminant  plumes  must  be  incorporated. 

10 


UNCLASSIFIED 


UNCLASSIFIED 


Most  of  the  processes  that  have  been  modeled  using  the  Lagrangian  approach  are 
straightforward  processes  such  as  linear  sorption  and  decay,  which  can  be  handled  in  a 
particle  analysis  by  adjusting  the  velocity  and  mass  of  the  particles.  More  complex  and 
coupled  interactions  can  be  approached  through  sequentially  solving  advective  and 
reactive  transport  equations.  There  are  examples  in  the  ground-water  literature  of  coupled 
transport  and  geochemical  models  (e.g.,  Gerla,  1992;  Fabriol,  1993;  Garcia-Delgado  and 
Koussis,  1997).  Streamline  methods  recently  used  for  modeling  a  variety  of  transport 
phenomenon  in  petroleum  reservoirs  provide  an  intuitive  approach  to  coupling  more 

complete  reactive  transport  models  with  ground-water  flow  models  to  provide  new  insight 
into  reactive  transport. 


e.  Streamline  Methods 


A  streamline  method  for  flow  in  porous  media  is  a  computational  technique  that 
approximates  the  transport  of  fluid  constituents  using  a  collection  of  one-dimensional 
mass  conservation  equations.  These  mass  conservation  equations  correspond  in  a  one-to- 
one  fashion  with  the  number  of  streamlines  in  the  model.  The  path  each  streamline  will 
follow  is  computed  in  the  same  way  as  the  tracing  of  pathlines  described  by  Pollock 
(1988)  and  discussed  earlier.  Therefore,  the  streamline  method  requires  the  solution  to  the 
ground-water  flow  equation  as  input.  Additionally,  the  formulation  of  a  streamline 
method  requires  that  the  appropriate  mass  conservation  equation  be  transformed  from  its 
three-dimensional  form  to  a  one-dimensional  form  along  a  streamline  This  transformation 
is  summarized  below  for  the  case  of  single-phase,  multicomponent,  incompressible  flow 
(after  Thiele  et  al„  1997).  In  this  case,  the  multidimensional  form  of  the  mass 
conservation  equation  can  be  expressed  as 

dC _ 

^  +u,'VFj-  0,  for  /=  1,....,  number  of  components,  (11) 


where  C,  is  the  concentration  of  component  i  and  F,  is  the  convective  flux  of  component  i. 
By  defining  a  time-of-flight  coordinate  along  a  streamline  as 


■  / 


<> 


s , 


where  s  is  the  streamline  direction,  the  divergence  operator  can  be 


(12) 

rewritten  according  to 
(13) 
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Substitution  in  Equation  1 1  transforms  the  conservation  equation  to  a  one¬ 
dimensional  formulation  for  each  component 

dC,  dFj  f 

=  0 ,  for  i  =  1,...,  number  of  components.  (14) 

A  one-dimensional  numerical  solver  is  used  to  solve  Equation  14.  This  method 
can  be  extended  to  any  equation  that  describes  the  physics  of  interest  in  one  dimension. 

A  streamline  method  for  transport  in  porous  media  is  better  suited  than  particle 
tracking  for  looking  at  the  complexities  of  reactive  transport.  Once  streamlines  are  traced 
through  a  three-dimensional  flow  field,  each  streamline  is  treated  as  a  one-dimensional 
system  along  which  mass  is  transported  by  a  model  that  appropriately  captures  the  physics 
of  flow.  Current  reactive  transport  models  that  embody  the  most  complete  and  rigorous 
representation  of  chemical  processes  are  generally  limited  to  one-  and  two-dimensional 
saturated  flow  systems.  Full  three-dimensional  reactive  transport  simulations  are  typically 
very  coarsely  resolved  (Johnson  et  al.,  1997).  For  example,  Fabriol  et  al.  (1993)  used 
coupled  reactive  transport  simulator  with  a  five-by-five  grid  to  predict  the  composition  of 
water  that  percolated  through  sandstone.  Clearly,  a  fully  coupled  reactive  transport 
simulator  could  not  be  used  to  solve  field-scale  contaminant  transport  model  unless  an 
extremely  coarse  grid  is  used.  The  trend  in  modeling  complex  transport  processes  has 
been  toward  decoupled  or  sequential  solutions,  which  are  conceptually  intuitive.  These 

solutions  involve  coupling  together  flow,  advective,  dispersive,  and  reactive  transport 
models. 

The  streamline  approach  has  already  been  used  to  model  numerous  single-  and 
two-phase  transport  problems  in  heterogeneous  petroleum  reservoirs  (Batycky  et  al., 
1997).  Most  of  the  mass  transport  in  particle  tracking  considers  the  movement  of  a  single 
component  through  the  flow  field.  If  many  components  are  to  be  considered,  the 
accounting  issues  become  overwhelming.  Since  the  streamline  method  moves 
compositions  instead  of  mass,  dealing  with  multicomponent  systems  is  computationally 
more  efficient.  In  addition,  particle  tracking  is  not  well  suited  to  looking  at  coupling 
between  the  components  in  the  system.  Most  pollutants  in  ground  water  contain  multiple 
components,  which  react  with  one  another  and  the  porous  media.  These  components  are 
coupled  together  through  processes  such  as  radioactive  decay  or  competition  for  sorption 
sites.  Equations  such  as  those  describing  the  migration  of  a  radionuclide  decay  chain  can 
easily  be  formulated  as  a  set  of  coupled  one-dimensional  equations  (Gureghian  and 


12 

UNCLASSIFIED 


UNCLASSIFIED 


Jansen,  1985).  This  set  of  equations  can  be  solved  along  each  streamline,  thereby 
capturing  the  evolution  of  the  solute  over  time. 

It  seems  clear  that  the  field  of  solute  transport  modeling  is  moving  toward  a  more 
complete  description  of  the  many  processes  that  influence  the  development  of  contaminant 
plumes.  This  process  necessarily  involves  looking  at  full  three-dimensional 
representations  of  the  flow  field  and  solute  mass.  But  the  most  complete  reactive  transport 
models  are  one-dimensional.  Since  streamlines  are  one-dimensional  pathways,  the 
reactive  transport  models  can  be  applied  to  each  streamline  in  the  system.  Because  the 
intersection  of  the  streamlines  with  an  underlying  three-dimensional  grid  is  also  known, 
these  one-dimensional  reactive  transport  solutions  can  be  placed  in  a  three-dimensional 
setting.  This  could  greatly  enhance  the  ability  to  model  field-scale  contaminant  transport 
problems. 

3.  Conclusions 

Streamline-based  methods  recently  developed  to  simulate  multiphase  transport  in 
petroleum  reservoirs  may  be  extended  to  model  contaminant  transport.  This  method 
overcomes  the  limitations  of  particle  tracking  techniques,  and  the  errors  and  poor 
performance  of  conventional  grid-based  methods.  Using  streamline  methods,  there  is  the 
potential  to  enhance  considerably  the  ability  to  model  complex  nonlinear  transport 
phenomena  in  heterogeneous  aquifers. 
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